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Magnetars are neutron stars with extremely high magnetic fields (= 


1014 Gauss) that exhibit various X-ray phenomena such as sporadic sub- 


second bursts, long-term persistent flux enhancements, and variable rates of 


rotation period change [1, 2]. In 2020, a fast radio burst (FRB), akin to cos- 


mological millisecond-duration radio bursts, was detected from the Galactic 


magnetar SGR 193542154 [3-5], confirming the long-suspected association 
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between some FRBs and magnetars. However, the mechanism for FRB gen- 
eration in magnetars remains unclear. Here we report the X-ray discovery of 
an unprecedented double glitch in SGR 1935+2154 within a time interval of 
approximately nine hours, bracketing an FRB that occurred on October 14, 
2022 [6, 7]. Each glitch involved a significant increase in the magnetar’s spin 
frequency, being among the largest abrupt changes in neutron star rotation 
[8-10] ever observed. Between the glitches, the magnetar exhibited a rapid 
spin-down phase, accompanied by a profound increase and subsequent decline 
in its persistent X-ray emission and burst rate. We postulate that a strong, 
ephemeral, magnetospheric wind [11] provides the torque that rapidly slows 
the star’s rotation. The trigger for the first glitch couples the star’s crust to 
its magnetosphere, enhances the various X-ray signals, and spawns the wind 
that alters magnetospheric conditions that might produce the FRB. 

SGR 1935+2154, the most active Galactic magnetar of the last decade 
with a spin period of 3.25 s and a dipole (polar) magnetic field strength of 
4 x 1014 G [12], has shown several major X-ray and radio outbursts [13, 14]. 
The most recent outburst episode of SGR 1935+2154 occurred in October 
2022 and lasted several days during which it emitted hundreds of short X-ray 
bursts [15-17] (see below). At Coordinated Universal Time (UTC) 19:21:47 
(topocentric time of the Canadian Hydrogen Intensity Mapping Experiment, 
hereafter CHIME) on October 14, 2022, this magnetar emitted an FRB with 
multiple radio peaks measured with CHIME and the Robert C. Byrd Green 
Bank Telescope (GBT) [6, 7]. 

Alerted prior to this activity [15], we initiated a series of Neutron Star Inte- 
rior Composition Explorer (NICER) with 75 ks exposure and Nuclear Spec- 
troscopic Telescope Array (NuSTAR) X-ray observations with 96 ks exposure 
starting from 17:32:40 UTC on October 12 until November 06, 2022 (Extended 
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Data Table 1). Hereafter, we use the elapsed time t since the FRB detection 
at barycentric corrected Modified Julian Date (MJD) 59866.80817034. Our 
high-cadence X-ray observations attained 67% temporal coverage bracketing 
the FRB, —17 < t < 10 hours, but did not cover the FRB itself due to Earth 
occultation. Nevertheless, this cadence enabled an unprecedented investigation 
of the spin-frequency evolution of the magnetar around the FRB. After remov- 
ing short bursts with timescales from milliseconds to around ten seconds, the 
stellar rotational frequency v can be revealed by pulsed thermal emission from 
the stellar surface using NICER in 2-8 keV and NuSTAR in 3-8 keV. Initially, 
the spin behavior exhibited a spin-down rate of ù = —1.7(6) x 10-11 Hz s7! 
(Figure la), with parentheses indicating the lo uncertainty on the last digit. 

A precise timing model was then derived by phase-coherent analysis, i.e., 
tracking the X-ray pulse times-of-arrival (TOAs) with a predictive model and 
optimizing its parameters [18]. In the phase evolution (Figure la and Extended 
Data Figure 1), a concave trend is observed in the range of —4.4 < t < 
4.4 hours, followed by a subsequent long-term convex quadratic trend. These 
data points are satisfactorily fitted by a model (see Methods) incorporating 
two spin-up glitches at tj, = —4.419°3 hours (MJD 59866.6310'}5) and tg2 = 


44105 hours (MJD 59866.99 + 0.02). With an uncertainty of approximately 


+30 minutes, these constitute the best-constrained glitch epochs ever observed 
in magnetars and pulsars at high energies [8, 19]. We attempted to fit the 
TOAs using a high-order polynomial but the resulting reduced x2 value was 
much higher than that of the best-fit two-glitch model (see Methods). 

At the first glitch, the frequency jumps by Av; = 3.0(3) x 107° Hz, 
corresponding to a fractional change of Av;/v = 1.0(1) x 1074. The spin- 
down rate ù also changed by Av, = —1.5(3) x 107° Hz s~!. The second 
glitch had a smaller Avy = 1.9(3) x 107° Hz (Av2/v = 6(1) x 1075) with 


Springer Nature 2021 TeX template 


SGR 1935 FRB and Glitch 5 


Av, = 1.5(3) x 107° Hz s71, thereby compensating for the excess recovery that 
occurred between the two glitches. Uncertainties are obtained through Markov- 
Chain Monte Carlo (MCMC) simulations (see Methods). Figure 2 compares 
the detected glitches with previously measured magnitudes in magnetars and 
pulsars [8, 10, 11, 20, 21], demonstrating that these two glitches are among the 
largest ever observed. The rotational energy increments associated with these 
two glitches are 3.9 x 1011 erg and 2.6 x 1041 erg, respectively (see Methods). 

Panels b and c of Figure 1 show the frequency evolution and its derivative 
using the two-glitch model. Additionally, we obtained the local spin ephemeris 
by allowing overlaps among consecutive TOAs and avoiding time intervals 
that cross the glitch epochs. The |/| between the two glitches increased by a 
factor of ~ 100 compared to the pre-glitch value, signaling a rapid spin-down 
bracketing the time of the FRB. 

To investigate the evolution of short burst occurrence rate (Figure 1d) 
in this episode, we applied the Bayesian-block technique to search for burst 
candidates using the photon events collected by NICER and NuSTAR (see 
Methods). The burst-active epoch lasted for at least two days, during which the 
rate remained at approximately one burst per minute. Subsequently, at t x —2 
hours, the rate peaked at more than four bursts per minute and subsequently 
dropped to below 0.05 bursts per minute in the next ten hours, being less than 
10% of the pre-glitch level at t42. This declining phase is contemporaneous with 
the occurrence of the FRB. Such interplay between the X-ray bursts and FRB 
is similar to the behavior seen during the previous outburst of SGR 1935+2154 
observed in 2020 [14, 22]. 

After excising these burst intervals, we extracted the X-ray photons and 
plotted the 64-s binned count rates of the persistent X-ray emission (Figure 


le). The NICER and NuSTAR count rates at the early outburst phase (—51 < 
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t < —2.5 hours) were approximately 1.0 and 0.5 count s~!. The persistent X- 
ray flux remained constant until a sudden increase at t = —2.5 hours, reaching 


1 


a NuSTAR count rate of over 50 count s within one minute. Subsequently, 


1 within ten hours. We 


the count rate declined to approximately 0.5 count s7 
also note that the FRB is phase-correlated with the 20-79 keV pulse peak, 
which exhibits a 0.5 phase offset relative to the soft X-ray pulse profile (see 
Methods). 

We also investigated the NICER and NuSTAR X-ray spectral behavior 
contemporaneous to the FRB (Figure 3). We divided the data into six time 
periods; (A) the early stage t < tz, (B) from t,; < t < —2.5 hours, the 
peak of the X-ray emission (C) —2.5 < t < —0.2 hours before the FRB, (D) 
0.5 < t < tg2, (E) decaying phase at tg2 < t < 10 hours (Figure la), and (F) 
the returning to quiescence at 11.5 days < t < 12.2 days (Extended Data Table 
3). The best-fit spectral parameters are summarized in Extended Data Table 4. 
The quiescent spectrum is well characterized by two components, as reported 
in other magnetars [23]; a blackbody from a hotspot on the stellar surface 


with a temperature of kT = 0.427100 keV and a radius R = a km when 


assuming a fiducial distance d = 6.6 kpc [24], and a non-thermal power-law 


component from the magnetosphere with a photon index of I = 1.2 + 0.1. 
At epoch A, the persistent X-ray flux increased by a factor of five com- 
pared to the quiescent flux level. The spectrum is modeled by the same 


two-component model but with a hotter blackbody temperature at kT = 


0.51+0.03 keV with R = 2.1+0.2 km and I = 1.1+0.1. At epoch B, no signif- 
icant flux change was observed. Subsequently, the flux increased by a factor of 
> 10 at 2.5 hours before the FRB (epoch C), with the power-law component 
dominating the X-ray spectrum (see detailed analysis in Methods). After the 


FRB (epoch D), the flux gradually decreased and, just after the second glitch 
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(epoch E), returned to a similar value as that of epoch A with a soft I of 1.6- 
1.8. At epoch F, the flux dropped to half of the level in epoch E, yet was still 
3 times higher than the quiescent flux. 

The total isotropic-equivalent radiative energy of the persistent emission is 
estimated as 5.0 x 104° (d/6.6 kpc)? erg between the two glitches. The burst 
emission at epochs B, C, and D has a total energy of 2.3 x 10*° (d/6.6 kpc)? erg 
when fitted with a power-law spectrum. This X-ray radiative output, 7.3 x 
104° (d/6.6 kpc)? erg, is one order of magnitude lower than the spin-down 
inferred energy but two orders of magnitude higher than the normal quiescent 
energy of < 1029 (d/6.6 kpc)? erg within the same time span (Extended Data 
Table 4). 

We outline one possible picture of the unique behavior detailed above, 
which may provide clues to the environment that permits or triggers the 
observed FRB, and radio bursts from SGR 1935+2154 in general. Owing to 
the absence of major changes in the count rate, pulse profile, and spectra 
around tji, the first glitch was probably triggered by a mechanism internal 
to the magnetar. Spin-up glitches in normal pulsars [8-10] and some magne- 
tars [19] are thought to be due to angular momentum transfer from superfluid 
neutrons in the rapidly rotating inner crust (and possibly core) to the other 
parts of the star [25], with the latter being slowed over long times by magneto- 
spheric torques. The rapid spin-down after the first glitch quickly regenerated 
the spin lag between the superfluid and the rest of the star, leading to the sec- 
ond glitch. Using angular momentum conservation arguments (see Methods), 
the short time interval of 8.8 hours between glitches suggests a superfluid com- 
ponent of several tens of percent. This implies that a large fraction of the core, 


along with most free neutrons in the inner crust, is in a superfluid state. 
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Long-term magnetic activity in the outer layers can provide heat or excite 
oscillations that can induce the unpinning of superfluid vortices in deeper 
layers [26, 27]. The motion of a large number of superfluid vortices could move 
superconducting flux tubes in the core. The latter would then change the 
surface layer magnetic field geometry, stress the crust, and heat its outer layers, 
thereby possibly producing the increased persistent emission and burst rate a 
few hours after the first glitch (see Methods). These stresses rupture the crust 
near the magnetic poles where the field is strongest. The associated heating 
to super-Eddington temperatures can lead to the expulsion of large amounts 
of ion-rich plasma (“volcanism” into the magnetosphere), forming an optically 
thick and collimated relativistic wind that torques the star. This ephemeral 
wind can naturally explain the rapid spin-down epoch between the glitches, 
provided the cumulative mass expelled is sufficient (see Methods). The wind 
twists field lines and increases the magnetic field energy density [28, 29] in the 
magnetosphere over the poles. As it abates, it untwists the field and clears the 
magnetosphere of current, perhaps seeding the decline of the X-ray flux and 
burst rate before the radio burst (see Figure 1d, le). 

A possible scenario for the FRB is the ignition of pair cascades by a mag- 
netar short burst as the magnetosphere evolves to a more charge-starved state 
[30]. Such phenomenology is generally consistent with the events around the 
April 2020 FRB, which occurred many hours into the waning phase of intense 
burst activity [14] and the post-anti-glitch radio burst in October 2020 [11]. 
Future high-cadence observations of SGR 1935+2154 and other magnetars in 
X-rays, like those presented here and in conjunction with radio monitoring, 


will help identify the conditions required for generating FRBs. 
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Fig. 1: Evolution of pulse timing, short X-ray bursts, and flux near 
the double-glitch epoch. The time zero epoch, marked by a vertical red 
dashed line, is defined as the occurrence time of the CHIME FRB. The dashed- 
dotted and dotted vertical lines denote the times of the first glitch and second 
glitch, respectively. Panel a shows the evolution of the pulse arrival phase. 
The orange line denotes the best-fit timing solution and the gray-filled area is 
the 1-o confidence interval obtained with MCMC simulations. Black and gray 
horizontal lines (labeled A through E) denote the time intervals of spectral 
analysis. Panel b describes the frequency evolution derived from the best- 
fit timing model (orange line) and the corresponding 1-ø confidence interval. 
The data points are best-fit frequencies obtained from 15 consecutive arrival 
phases with a moving window of one data point. Panel c is the ù evolution 
obtained using the same procedure as in panel b. Panel d shows the evolution 
of the burst occurrence rate, where blue circles denote the values obtained 
with NICER and green squares denote those obtained with NuSTAR. Panel e 
shows the count rate of the persistent emission after removing the bursts. 
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Fig. 2: Pulse frequency discontinuities for known pulsar glitches. 
Glitch size distribution of previously reported rotation-powered pulsars (cir- 
cles) and magnetars (squares) [8, 10] in the phase space of changes of stellar 
spin frequencies (ordinate) and their derivatives (abscissa) between glitches. 
The data are shown as absolute values. The glitches of rotation-powered pul- 
sars are bimodally distributed, where glitches of Vela-like pulsars (blue circles 
filled with pink) mainly lie on the top of the figure, and glitches of other pul- 
sars form a linear trend with scatter [8]. The two magnetar glitches observed 
in this work are indicated by red diamonds. They have extremely large changes 
in |Av| compared to all pulsars, and |Av| comparable to those of Vela-like 
pulsars. These two glitches lie on the extension of the linear trend of pulsar 
glitches. 
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Fig. 3: X-ray spectral variations near the double-glitch epoch. X-ray 
spectral evolution of the persistent emission (red) of SGR 1935+2154 after 
removing short bursts, obtained with NICER and NuSTAR in time epochs A to 
F (see panel a of Figure 1). NICER spectra are shown in crosses, and NuSTAR 
FPMA and FPMB spectra are shown in circles and triangles, respectively. As 
a reference, we compare each of these active-state spectra with one quiescent 
spectrum (gray) obtained on October 16, 2020, two years before the present 
FRB. The spectra of epochs A, B, D, E, F, and the quiescent state are fit 
using a blackbody plus a power law with absorption in the interstellar medium. 
The spectra in epoch C are approximated with a power law. The spectral 
components, including photo-absorption, are plotted using dashed (blackbody) 
and dotted (power law) lines. 


Springer Nature 2021 TeX template 


16 SGR 1935 FRB and Glitch 


Methods 


NICER data reduction. NICER is a non-imaging X-ray observatory with 


2 at 1.5 keV and an excellent time reso- 


a large collecting area of 1900 cm 
lution of < 300 nanoseconds aboard the International Space Station [31]. It 
consists of 56 co-aligned X-ray concentrating optics and associated silicon- 
drift detectors, of which 52 are in operation. We reprocessed the NICER data 
using NICERDAS version 9, part of the HEASOFT version 6.30.1 package 
and the calibration database version 20210720. The unfiltered data were cali- 
brated and screened using the nicerl2 pipeline with default settings, resulting 
in clean scientific event lists. We further removed a few background flares 
near GTI boundaries. We corrected the arrival time of each photon to the 
barycenter of the solar system using the barycorr command with the JPL 
solar-system ephemeris DE405. Photons with energies below 2 keV, due to 
significant absorption, and above 8 keV, due to reduced NICER sensitivity, 
were excluded from the analysis. A summary of the NICER data used in this 
research can be found in Extended Data Table 1. 

NuSTAR data reduction. NuSTAR consists of two identical focal plane 
modules (FPM), FPMA and FPMB. Each module is equipped with conical 
Wolter-I optics that achieve a point spread function with a half-power diam- 
eter of around 60 arcseconds. Each FPM houses a detector with an energy 
range of 3-79 keV, a field of view of 12’ x 12’, an energy resolution of 0.4 
keV at 10 keV, and a timing resolution of 65 microseconds[32]. The data 
were reduced using NUSTARDAS version 2.1.2 with the calibration database 
version 20211020. We performed data calibration and screening using the 
nupipeline script. Images, event lists, light curves, and spectra were extracted 
using nuproducts. The source photons were extracted from a circular region 


with a radius of 0.85 arcmin centered on the position of SGR 1935+2154. This 
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guarantees that around 85% of the source X-ray photons were encircled. We 
extracted the background spectra from a circular source-free region within a 
radius of 2.5 arcmin centered on R.A. 19:34:46.6029 and Decl. +21:46:55.710 
in the J2000 reference frame. 

Burst identification: We used the Bayesian block technique!, which is 
designed to detect local variability by dividing data into variable-sized time 
blocks, to search for significant flux changes [33]. The search energy ranges for 
the NICER and NuSTAR events were 2-8 keV and 3-79 keV, respectively. The 
prior was computed with a false alarm probability of po = 0.01. We first iden- 
tified blocks with a duration shorter than 4 seconds, and then merged nearby 
blocks as bright and long bursts are typically identified as several consecutive 
blocks. Next, we assessed the likelihood of the total counts in each block being 
a Poisson random fluctuation around the count rate present approximately 
1 second before and after the block. We considered all bursts with a proba- 
bility lower than 3 x 1077 (corresponding to a detection significance of 50) to 
be confirmed bursts and marked the rest as burst candidates. We removed all 
the confirmed bursts and burst candidates from the analysis of the persistent 
count rate and rotation frequency. In addition, NuSTAR captured a flare with 
multiple peaks that occurred approximately two hours before the time of the 
FRB emission. The persistent emission during the peak of the flare increased 
by two orders of magnitude and shows multiple peaks, making it challenging to 
identify individual bursts. We do not attempt to identify the numerous small 
peaks within the 200-second long flare. We conservatively count it as a single 
burst and remove it from further timing analysis since the persistent emission 
between peaks may be heavily contaminated. Therefore, the burst rate at the 
flux peak in Figure 1d may be underestimated. A total of 368 short burst can- 


didates were identified in the NICER data set and 720 in the NuSTAR data 


https: //docs.astropy.org/en/stable/api/astropy.stats. bayesian_blocks.html 
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set. Before the FRB, SGR 1935+2154 exhibited a high burst rate of > 1 short 
bursts per minute, which then dropped to < 0.1 bursts per minute within a 
time scale of ~5 hours. 

Timing analysis and glitch model: We first collected X-ray photons 
detected with NICER before the FRB and performed a two-dimensional Z?- 
test with two harmonics (n = 2) [34]. As the time span increased, we stopped 
the Z2 search when the spin-down rate could be significantly determined, 
which was v = 0.3075277(1) Hz and ù — 2.0(6) x 1071 Hz s~! at a time zero 
epoch of MJD 59865.678 (TDB). This corresponds to a spin period derivative 
of P = 2.1(6) x 107!° ss~!, an order of magnitude higher than the long- 
term value of P = 1.43(1) x 1071! ss~! [35, 36]. This result suggests that the 
spin behavior of SGR 1935+2154 significantly changed after the initial activity 
detected on October 10 [15]. 

To ensure high accuracy in tracing the TOAs, we divided the NICER obser- 
vations into segments of integer orbits. Each segment was chosen to contain 
a sufficient number of photons to give a small uncertainty of og < 0.05 for 
each TOA. Near the peak of the flare, a NICER orbit was further divided into 
three segments owing to a large number of collected photons. The TOAs were 
determined using the unbinned maximum-likelihood method [18, 37]. A pulse 
profile with the highest signal-to-noise ratio was selected to construct a non- 
parametric template by smoothing it with a Gaussian filter and interpolating 
it to a bin size of 1074 cycles. This same step size was used to compute the 
phase shift A% of each TOA. The same analysis procedure was also performed 


on the NuSTAR data in the 3-8 keV energy range. 
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Then, we traced the phase evolution of the TOAs based on the aforemen- 


tioned ephemeris (see Extended Data Figure 1): 


A(t) = v(t — to) + = D(t — to)? 4 s(t to) +... (1) 


where those terms higher than the first derivative term (ù) are ignored for 
the first trial. It was clear that the initial timing solution failed to accurately 
predict the subsequent evolution. The phases of the TOAs showed a drop 
of approximately 0.2 cycles over a timescale of S 0.5 days, starting about 
4.4 hours prior to the CHIME FRB detection. Interestingly, the evolution 
formed a concave quadratic curve, indicating a fast spin-down over half a 
day. After approximately 4.4 hours following the CHIME FRB, the evolution 
showed a long-term convex quadratic curve, implying a much slower spin-down 


thereafter. Consequently, we fit the TOAs with a two-glitch model: 


P(t) = polt) = v(t — to) + 50(t—to)? for t< tg 

o(t) = p1 (t) = polt) + Avı (t — tA) + FAM (t—tyi)? for ta <t < tye 

Plt) = p1 (t) + Ano(t — tg2) + 4Ai2(t — tg2)? for t > tye 

(2) 

where t,; and tg2 are glitch times. The sudden changes in frequency and spin- 
down rate are denoted as Avı and Aù; for the first glitch, and Avg and Aro 
for the second glitch, respectively. The best-fit parameters for this model are 
listed in Extended Data Table 2. The best-fit timing model is shown as the 
orange curve in Figure 1 and Extended Data Figure 1. This model yields a fit 
statistic x? value of 83.2 with 64 degrees of freedom, and the residual did not 


show significant long-term trends. 
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The best-fit parameters and their corresponding 1-o uncertainties were 
obtained through an MCMC simulation. We sampled the full parameter space 
of the two-glitch model using the encee MCMC sampler[38] and assuming flat 
prior probability densities for all parameters. The MCMC was run with 32 
walkers for 10° steps, and the first 101 samplings were discarded. The posterior 
probability density distributions are displayed in Extended Data Figure 2, and 
the values and 1-o uncertainties are summarized in Extended Data Table 2. 
The 1-c intervals of the best-fit model are plotted as the gray shaded area in 
Figure 1 and Extended Data Figure 1. 

In comparison to the pre-glitch epochs, the post-glitch ephemeris is more 
precisely defined owing to its long time span. To verify our timing solution, 
we cross-checked it by calculating the phases of the TOAs based on the post- 
glitch second-order ephemeris with a reference epoch of to =MJD 59868.87136, 
v = 0.3075264 Hz, and ù = —4.9x 10712 Hzs~!. A zoom-in view of the phase 
evolution near the FRB is presented in Extended Data Figure 3. The parame- 
ters and the x? statistic of the best-fit model are in full agreement with those 
obtained previously. Furthermore, we tried to fit the phase evolution with an 
8t_order polynomial, including terms up to dv/dt? in Equation 1. Although 
this high-order polynomial has the same number of free parameters as the 
two-glitch model, it results in a much higher x2 = 378.8 with 64 degrees of free- 
dom. The residuals (as shown in panel c of Extended Data Figure 3) exhibit 
a noticeable long-term trend, implying that the two-glitch model provides a 
much better fit. 

The two glitches observed in this outburst are among the largest of all 
glitches determined in pulsars and magnetars, as demonstrated in Figure 2. We 


plotted the values of Av versus Av and compared them to previously recorded 
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glitches in pulsars and magnetars obtained from the literature and the Jodrell 
Bank Glitch Catalogue [8, 11, 20, 21]. 

The rotational energy change involved in a glitch is Eglitch © An? IvAv, 
where I (~ 1045 g cm?) is the stellar moment of inertia. Thus the energies 
associated with these two glitches are 3.9 x 1011 erg and 2.6 x 1041 erg. The 


change in rotational energy between the glitches is estimated to be 


AE vot = ~ Exot dt œ —6.5 x 1011 erg, (3) 
tgi 
where His = 4n?Ivv is the inferred rotational energy loss rate during the 
period of enhanced torque. This value of AZ,,; is approximately equal to the 
total energy associated with the two glitches. 
Energy-Resolved Pulse Profile: We investigated the energy dependence of 
the pulse profile to gain insights into the emission mechanisms. To this end, 
we obtained the energy-resolved pulse profiles based on the best-fit two-glitch 
model by dividing the NICER data into two energy bands, 2-5 keV and 5- 
10 keV, and the NuSTAR data into four energy bands, 3-5 keV, 5-10 keV, 
10-20 keV, and 20-79 keV. The results are presented in Extended Data Figure 
4. The soft X-ray peak is evident at a phase of ¢ ~ 0, while the hard X-ray 
peak is found at % ~ 0.5. The soft X-ray peak dominates the profile below 10 
keV but vanishes at energies above 20 keV. On the other hand, the hard X-ray 
peak is present above 5 keV and becomes the dominant feature above 10 keV. 
We also calculated the spin phase of the FRB. The barycentric cor- 
rected time of the CHIME FRB[6] was calculated to be MJD 59866.80817034 
using pintbary, a command line tool of the pulsar timing software package 
pint|39]2. Similarly, we calculated the barycentric correct time of the brightest 


of the GPT FRB[7] as MJD 59866.80817889. These times were extrapolated 


*http://ascl.net/1902.007 
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to infinite frequency based on a dispersion measure of 332.8 [6]. The phase of 
the FRB was found to be approximately 0.38 based on CHIME detection or 
0.5 based on GBT result, which is close to the peak of the 20-79 keV band 
(see Extended Data Figure 4). 

Verification of the spin-down rate between glitches: From the analy- 
sis of the energy-resolved pulse profile, it is evident that the pulsed emission 
has multiple origins. The soft X-ray band may be dominated by thermal emis- 
sion from the neutron star’s surface, whereas the hard X-ray emission would 
be originating from the magnetosphere. This suggests that the timing solution 
based on photons below 8 keV could be affected by the pulse profile change 
due to the motion of hotspots on the surface during the outburst. To investi- 
gate this, we attempted to trace the evolution of the phase of the 20-79 keV 
profile, but the large uncertainty due to the limited number of photons col- 
lected (approximately 2200 counts between the two glitches) made it difficult 
to trace the evolution in detail. 

As an independent verification, we performed an energy-resolved two- 
dimensional Z3-test search for v and £ between the two glitches. Using 
NuSTAR events collected between the two glitches, we divided the events into 
four energy bands as defined in the energy-resolved pulse profile analysis. The 
number of photons in each energy band was 7991 in the 3-5 keV range, 14478 
in the 5-10 keV range, 6484 in the 10-20 keV range, and 2275 in the 20-79 
keV range. The results of the two-dimensional Z2-test for each energy band, 
as well as the broadband (3-79 keV) events, are presented in Extended Data 
Figure 5. The result obtained from the broadband events suggests a spin-down 
rate of ù = —1.7(4) x 107° Hz s~!, which is in full agreement with the value 
obtained from the TOA analysis [—1.5(3) x 107° Hz s~']. The spin-down rate 


in each energy band was determined to be —1.9(9) x 107° Hz s~! for 3-5 keV, 
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—1.5(5) x 107° Hz s~! for 5-10 keV, —1.8(7) x 107° Hz s~! for 10-20 keV, and 
—1.3(9) x 107° Hz s~! for 20-79 keV. These results point to a large spin-down 
rate that is consistent with that obtained from the TOA analysis. 

Moreover, we obtained the soft X-ray profile 2-8 keV NICER data and 10- 

79 keV NuSTAR data in three epochs: t < tg1, tgi < t < tg2, and t > tga. Those 
X-ray photons collected one hour close to the glitch epochs are not included 
to avoid possible coverage of the glitch time. No significant phase shift can 
be observed. These suggest the recovered phase evolution is secure and not a 
misinterpretation of hotspot migration. 
Spectral Analysis: We extracted NICER and NuSTAR (FPMA and 
FPMB) persistent spectra with burst-removed and barycentric corrected GTIs 
(Extended Data Table 3). NICER spectra were extracted with xselect version 
2.5a following the standard procedure, where the background of NICER data 
was produced with NICER background estimator tool nibackgen3C50 version 
v7b [40]. The persistent spectra of epochs A, D, and E are extracted from the 
time intervals where NICER and NuSTAR observed the target simultaneously, 
whereas those of epochs B, C, and F do not strictly rely on time intervals 
with simultaneous observations. We also extracted the burst spectra of NICER 
and NuSTAR (FPMA and FPMB) by accumulating all the detected burst at 
epochs B, C, and D using the barycentric corrected burst GTIs. 

We used the XSPEC software version 12.12.1 with the Chi-square statistic in 
the following analysis. We fitted the persistent spectra at each epoch (Figure 
3) and the burst accumulated spectrum (Extended Data Figure 6), where the 
information of observation and the best fit spectral parameters are shown in 
Extended Data Table 4 and 3, respectively. The absorbed blackbody (BB) 
and power law (PL) model, (tbabs* (bbodyradt+pegpwrlw) ), was used for the 


persistent spectra except for epoch C. The solar abundance value was used 
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to calculate the absorption [41]. An absorbed PL model (tbabs*pegpwrlw) 
was used for epoch C since the BB component is negligible compared to the 
PL. The accumulated burst spectrum was fitted with an absorbed BB plus 
a PL with a high energy exponential roll-off (tbabs* (bbodyrad*cutoffpl1) ), 
hereafter CPL. We added a constant term (~ 3) in the spectral model to cross- 
calibrate the normalization between NICER and NuSTAR because these two 
instruments cover emissions from different bursts. Best-fit spectral parameters, 
including the hydrogen column density Ny, BB temperature kT, its radius R, 
and the photon index I of the PL, are tabulated in Extended Data Table 4. 
These parameters were set free except for the spectra of epoch B, where the 
BB parameters are frozen to those of epoch A since no NICER observation 
was available and NuSTAR spectral shape was similar to that of epoch A than 
that of epoch C. 

The absorbed 2-70 keV X-ray flux of persistent spectra was calculated 
with the flux command of XSPEC and the NICER and NuSTAR spectra were 
used for low (2-10 keV) and high (10-70 keV) energy bands, respectively. 
The absorption-corrected flux of persistent spectra was calculated from the 
normalization factors of the BB and PL models. The absorbed and absorption 
corrected flux of burst spectra was calculated from normalization factors of 
the BB and PL models with only NuSTAR spectra because of the difficulty of 
estimating the effect of the constant value. 

The hydrogen column density Ny was determined to be (2.1-3.4) x 
1022 cm~? from the persistent spectra. Comparing the quiescent spectra, the 
BB temperature kT increased to 0.51-0.68 keV around the FRB (epoch A- 
E) and returned to the quiescent value 12 days after the FRB (epoch F). On 
the other hand, the BB radius R was similar to the quiescent value. Around 


the peak of the count rate in X-ray bands (Epoch C), the power-law spectrum 
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became softer (from IT ~ 1.2 to 1.6) and extended from the hard X-rays to 
lower energies, making it difficult to measure the BB parameters. In contrast, 
burst spectra have a high hydrogen column density (4.210:2 x 1022 cm7?) and 


a soft photon index (1.40 + 0.02). 


We also compared the above fitting results around the glitch (epochs A — 
E) with those considering SCORPEON background model (obtained with the 
tool nicerl3-spec of NICERDAS v10). A few fitting parameters vary from 
the values obtained with the 3C50 background model, but the difference is at 
most around 20% (see Extended Data Table 4). The difference could be larger 
than their lo statistical uncertainty and could be inherent in modeling the 
background. However, they did not change our order-of-magnitude estimate of 
the total energy and the spectral evolution. 
Theoretical Interpretation: The main deliverables of this unprecedented 
observational campaign are the two spin-up glitches that bracket a strong spin- 
down epoch that persists for several hours, the enhanced soft and hard X-ray 
activity evident mostly between the glitches, and the emission of the FRB 
during this activity interval. All these arise during an active epoch of X-ray 
bursts. The ephemeral wind picture presented here to address these features 
is derived from that presented in [11] for the spin-down glitch observed from 
SGR 1935+2154 in 2020 during a period where source monitoring was much 
sparser than for this event. 
Constraining the Superfluid Component: First, the focus is on the timing 
results. Spin-up glitches impart abrupt re-adjustments to the rotational equi- 
librium of the entire star [25, 42]. Large glitch sizes typically require moments 
of inertia Isf of the superfluid neutron component in the inner crust of around 


a few percent of that (T) of the entire star [43—46]. This ratio can be estimated 
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using angular momentum conservation arguments [43, 47], yielding 


Ma (a) 


I ` [ATs 

where AT, is the time interval between glitches and is ~ 8.8 hours in the case 
here. This estimate suggests a superfluid component comprising several tens 
of percent of the entire star. 

The Inter-Glitch Torque driven by a Strong Wind: The rapid spin-down 
episode following the first glitch leads to a net increment Av ~ —4.7 x 107° Hz 
in the hours prior to the second glitch. This evolution could be attributed to 
an external torque on the star and likely linked to the magnetosphere. The 
October 2020 anti-glitch (Av = —1.8x10 Š Hz) detailed in [11] was interpreted 
as due to a strong, ephemeral wind lasting around 10 hours or less in its 
peak phase. This wind loads portions of the magnetosphere that are tied to 
the active surface region with plasma, altering the magnetic configuration and 
opening field lines near the pole [48]. Thus the wind introduces a monopolar 
component to the field geometry that renders the polar fields more radial, 
approaching that of the morphology and strength of a split monopole [28, 29, 
49] in the limit of extremely large winds and strong twists. 

Such a picture can be also envisaged for the SGR 1935+2154 timing data 
for October 2022 presented here. The expected crustal stresses associated with 
the first glitch, amplified by the superstrong magnetic field, could precipitate 
plastic flows or rupturing, and localized heating and mechanical transport that 
could eject plasma from the surface near the magnetic pole into the magneto- 
sphere. There, the locally zenith-directed fields permit rapid energy and mass 
transport vertically, and the thermal radiation pressure is above the Edding- 


ton limit, driving baryon-loaded plasma into the magnetosphere as a wind. 
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Only mass transfer on open field lines can shed angular momentum to infinity 
and slows down the crust’s rotation. 

The wind required to drive the high spin-down episode between the glitches 
is highly optically thick [11], yet quite confined in its magnetic colatitudes in 
the inner magnetosphere. The polar cap colatitude of the last open field line 
footpoint at the surface for a dipolar configuration is cap © s/2mRns/ Pe N 
0.46° for a rotational period of P = 3.25 sec. Thus, even with a modest expan- 
sion of the polar cap due to twists at the surface, the wind collimation at 
low altitudes is sufficient that it does not significantly obscure the surface soft 
X-ray emission: the inferred blackbody emission zone radii of Rpp ~ 2km in 
Extended Date Table 4 suggest that the soft X-ray hot spot subtends angles 
bbb > Rop/(RasV7™) ~ 7°. This is for a source distance of 6.6 kpc, which is not 
well constrained, thereby introducing some uncertainty into this Op, estimate. 
Near the surface, the wind is far less dense than the plasma that radiates the 
much more impulsive magnetar giant flares, and its motion is much less rela- 
tivistic [11], thereby not being subject to strong Doppler boosting. Accordingly, 
its effective temperature lies at or below the NICER band [11], indicating that 
the wind’s periphery may contribute to the soft X-ray signal, but not over- 
whelm it with hard X rays. Yet such a wind does not occult the persistent 
hard X-ray signals at relatively low altitudes on closed field lines (see below). 

The analysis of ephemeral mass loss delivered in [11] that leads to the con- 
clusion of high optical depth in the wind needs to be adapted for the high |Av| 
observed here between the glitches. Using total angular momentum conserva- 
tion, the time-integrated mass loss ôm for a star of mass M and radius Ry, 
can be expressed via the ratio óm/M ~ 2/5 (Rns/Req)” |Av|/v for a spin down 
increment Av in v; see Eq. (3) in [11] and associated discussion for details. 


Here Req is the radius of equipartition between the energy densities of the 
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flowing plasma and the magnetic field. For the dipole field geometry focused 
on in [11], one quickly infers Reg ~ 5Rns and a mass loss ratio ~ 2.4 x 107°. 
This high value is not a realistic estimate because the plasma loading of the 
polar magnetosphere is sufficient to enhance the local field strengths [29] sub- 
stantially when approaching a split monopole configuration, wherein the field 
strength scales roughly as 1/R with altitude R. This twisted field enhance- 
ment moves the equipartition radius to around Reg ~ 20 — 50Rps, leading to 
the estimate of a fraction óm/M ~ 2.4 x 1078 — 1.5 x 1077 of the stellar mass 
expelled during a nominal 8.8-hour wind phase. This constitutes a mass loss 
roughly on the scale of that expected for giant flares. To provide a more pre- 
cise estimate requires a detailed study of the dynamics of radiation-coupled 
plasma winds and angular momentum expulsion in magnetar magnetospheres; 
this will be the subject of a future study. 

This wind can persist as long as the heat/mechanical energy source near 
the crustal active zone is maintained, and then abates as the surface regions 
“heal,” and the polar magnetosphere untwists. Observationally this “volcan- 
ism” timescale is <S 9 hours, and this provides an important constraint for 
theoretical models for how sub-surface conditions associated with magnetar 
glitches are generated. The cumulative spin-down during this epoch was suf- 
ficient to return the rotational configuration in ~ 9 hours to a critical state; 
this may be why a second glitch of similar magnitude ensued. The wind’s 
abatement suggests at least partial crustal stress relaxation or healing. Note 
that mass that is shed on peripheral closed field lines (spillover “wind” ) can 
potentially influence emission from the closed portions of the atmosphere and 


magnetosphere. 
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Enhanced X-ray Activity between the Glitches: Two hours after the 
first glitch, there is a sharp rise in the burst rate and the persistent hard X- 
ray emission flux from SGR 1935+2154 (see Figures 1 and 3). Assuming that 
these radiative changes are associated with the glitch, the time delay can be 
used to constrain the crustal location of heating within the magnetar, via the 
thermal timescale Ti, ~ (C/«)2z7, where C is heat capacity, « is thermal con- 
ductivity, and z is depth below the surface. For typical C and « of strongly 
magnetic neutron stars (see, e.g., [50]), a two-hour thermal timescale occurs 
at a depth z ~ 10 m and density of 107 g cm~? for polar locales where effi- 
cient heat transport occurs along B approximately in the radial direction. 
This timescale likely exceeds that for the onset of the wind, since this out- 
flow emanates from very near the surface; yet it may indicate when the wind 
enters its most powerful phase. The temporal correlation of the burst rate and 
persistent emission enhancement and decline following the first glitch suggests 
they are driven by the same sub-surface energy deposition, with concurrent 
evolution and relaxation of the crust, atmosphere, and magnetosphere. 
Concerning the spectroscopy, there may be a causal coupling between 
the observed heating of the surface soft X-rays in the NICER data and the 
enhanced persistent signal seen above 10 keV by NuSTAR between the glitches. 
The leading paradigm for the generation of the persistent hard X-ray emis- 
sion above 10 keV in magnetars is via resonant inverse Compton scattering 
(RICS) of soft X-ray photons emanating from the stellar surface by relativis- 
tic electrons accelerated in the magnetosphere [51, 52]. The scattering cross 
section is resonant at the cyclotron frequency. The electrons are accelerated 
by magnetospheric electric fields, coherent or non-stationary, and given that 
the NuSTAR signal is enhanced and not abated between the glitches, we can 


presume that the collimated wind does not impair lepton acceleration in the 
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pertinent electric field zone. This region is likely on closed field lines as most of 
the RICS signal, which is highly anisotropic and beamed, originates in quasi- 
equatorial locations where the magnetic field is low enough for the scattering 
to access the cyclotron resonance. Moreover, for these locations, an observer 
who samples the Doppler-boosted RICS emission along field tangents [53] will 
naturally see a hard X-ray pulse peak at a different phase from that for the 
soft X-ray surface signal, as is evident in Extended Data Figure 4. The inten- 
sity of the RICS emission is roughly proportional to the number density of 
the surface thermal photons [54], which for a blackbody scales as T’. Accord- 
ingly, the roughly 60% increase in T between the quiescent epoch and interval 
D as listed in Extended Data Table 4 should naively yield an enhancement 
in the RICS signal by a factor of ~ 4. In addition, the temperature increase 
might significantly extend the volumes of the inner magnetosphere for which 
the cyclotron resonance that controls the RICS emissivity is accessible, thereby 
further increasing the RICS flux. Combined, these influences would account 
for about half or more of the enhancement in the NuSTAR spectrum evident 
in Figure 3, the remaining portion perhaps due to increases in the relativistic 
electron density associated with a spillover “wind.” 

The Fast Radio Burst: If the electric fields that accelerate the electrons 
that generate the RICS signal are left unmodified by the polar wind, then their 
persistence may continue pair creation that contributes to an enhanced RICS 
signal that can then further lead to enhanced pair cascading. The absence or 
presence of plasma where a burst occurs is likely connected to the generation 
of an FRB [30, 55]. Thus, the timing of the FRB during the epoch of amplified 
and abating hard X-ray emission between the glitches may in fact be a causal 
association. As a putative spillover wind abated and the system relaxed after 


the local persistent flux peak 2-3 hours after the first glitch (time interval 
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C), this region became increasingly charge-starved to errant bursts during the 
span where the highest torque, and likely also stresses, were inflicted on the 
crust. This phenomenology is reminiscent of the April 2020 FRB, which also 
occurred in the waning phase of the burst storm on its track to relative source 
dormancy. Furthermore, it is noteworthy that the phase-connected arrival time 
of the FRB is roughly coincident with the peak of the hard X-ray signal in 
Extended Data Figure 4, and offset in phase with the peak of the soft X-rays. 
Thus it appears probable that in this case both the FRB and the persistent 
hard X-rays both originate from quasi-equatorial locations, though likely not 
at the same altitudes. 

The October 2022/October 2020 Glitch Comparison: It is of interest 
to compare this October 2022 event with the anti-glitch in SGR 1935+2154 
of the October 2020 epoch, which was not afforded the same high-quality 
observational cadence [11]. The 2020 anti-glitch was inferred using timing data 
spanning over 50 days, and arose during a data gap of approximately 3 days 
duration. An FRB-like burst was detected at about 3 days after the putative 
time of the anti-glitch, and pulsed radio emission ensued thereafter, both dur- 
ing a period when the X-ray timing was re-established and at the post-glitch 
ù [11]. No X-ray enhancement before or after the anti-glitch was seen. For 
the October 2022 events studied here, if instead, the source had been subject 
to a similar observational data gap of duration 1-2 days that bracketed the 
two glitches, none of the enhancement in persistent X-ray emission, the FRB 
and the glitches themselves would have been detected. In fact, the subsequent 
relatively sparse timing data might actually accommodate an anti-glitch solu- 
tion for the October 2022 event, albeit likely of different magnitude in Av/v. 
Clearly, the two events are different, particularly in terms of the eventual 


pulsed radio emission seen in 2020. Yet they do present similarities, principally 
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rapid spin-down episodes, and this motivates the position that an ephemeral 
wind can drive the strong torque on the star for several hours between the 
two glitches. Moreover, the events do pose the question of what sub-surface 
physics dictates a glitch recurrence timescale of around Tyee ~ 2 years in 
SGR 1935+2154. Glitch models of pulsars predict a range of Trec times and 
distributions [42]. Comparing these two events clearly indicates that the excel- 
lent data cadence we present here for October 2022 highlights the rich nature 


of information afforded by such temporally dense observations. 
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Extended Data Table 1: Data sets used in this work. 


Observatory ObsID Start Date Exposure (s) 
NICER 3020560159 2020-10-16T 21:40:00 1349 
NICER 5020560106 2022-10-12T17:40:27 1260 
NICER 5020560107  2022-10-13T00:03:40 14913 
NICER 5576010101 2022-10-13T16:44:01 3302 
NICER 5576010102 2022-10-14T03:34:37 10332 
NICER 5576010103 2022-10-14T 23:58:20 7757 
NICER 5576010104 2022-10-16T02:16:40 9990 
NICER 5576010105 2022-10-17T01:14:37 2633 
NICER 5576010106 2022-10-18T03:54:40 2472 
NICER 5576010107 2022-10-19 T04:20:01 2845 
NICER 5576010108 2022-10-20T06:58:40 1693 
NICER 5576010109 2022-10-21T00:00:00 3091 
NICER 5576010110 2022-10-22'T02:02:12 2034 
NICER 5020560108 2022-10-23 T21:44:37 260 
NICER 5020560109 2022-10-25T06:17:53 678 
NICER 5576010111 2022-10-26 T08:20:34 6592 
NICER 5576010112 2022-11-03 13:22:00 1297 
NICER 5576010113 2022-11-06 T08:16:30 1094 
NICER 5020560110 2022-11-05T09:05:40 1194 
NICER 5020560111 2022-11-06T00:42:00 729 

NuSTAR 90602332004  2020-10-16T21:16:09 18585 
NuSTAR 80802317002  2022-10-14T02:12:57 49624 
NuSTAR 80802317004  2022-10-26T08:35:11 46266 


Extended Data Table 2: Best fit spin parameters between 2022 October 12 
and November 06. Timing information for the glitches at times tj; and to, 
inter-glitch and post-glitch frequency derivative, and frequency shift between 
two glitches are also included. 


Parameter 


Value 


Epoch (MJD) 
u (Hz) 


59865.6782217 
0.3075277(6) 


ù (Hz s~') —1.7(6) x 1071! 
tgi (MJD) 59866.637 0-0; 
Av; (Hz) 3.0(3) x 1075 
Av (Hzs~') —1.5(3) x 1079 
pz (Hz s-1) —1.5(3) x 1079 
tg2 (MJD) 59866.99(2) 
Avy (Hz) 1.9(3) x 1075 
Avy (Hz s71) 1.5(3) x 107° 
pÍ (Hz s7!) —4.9(1) x 10712 
Av}, (Hz) 3.9(4) x 1075 
x? /dof 83.2/64 


* inter-glitch frequency derivative 
T post-glitch frequency derivative 
t frequency shift between two glitches 
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Extended Data Figure 1: Phase evolution of 'TOAs between 2022 October 
12 and November 06. The orange line is the best-fit two-glitch model, where 
the red vertical line denotes the time of the CHIME FRB detection. The times 
of two glitches are shown as the vertical dashed-dotted line (first glitch) and 
the dotted line (second glitch). The gray box denotes the zoom-in region shown 
in Figure 1. The residual is shown in the lower panel. 
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Extended Data Figure 2: Posterior probability density distributions of the 
emcee sampler The total simulation steps is 100,000. Two spin-up glitches are 
needed, where the first one is accompanied by an increase in the spin-down 
rate, while the second glitch is accompanied by a decrease in the spin-down 
rate. The change in the spin-down rate of these two glitches is mostly canceled 
out. In the 1-D histograms, the dashed lines represent the best-fit value along 
with its lo standard deviation. 
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Extended Data Figure 3: Phase evolution of TOAs determined with NICER 
and NuSTAR based on the post-glitch ephemeris. The symbols used here is the 
same as those in Extended Data Figure 1. The red dashed curve is the best-fit 
polynomial ephemeris up to the 8th order time derivative of v. This model has 
exactly the same number of free parameters as that of the two-glitch model. 
The residual of the two-glitch model is shown in panel (b), while that of the 
8th_order polynomial is shown in panel (c). 
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Extended Data Figure 4: Energy-resolved pulse profile of SGR 1935+2154.. 
Pulse profiles obtained with NICER are divided into 2-5 keV, 5-10 keV bands, 
where NuSTAR data were divided into 3-5 keV, 5-10 keV, 10-20 keV, and 
20-79 keV bands. The red dashed line and dotted line indicate the phase of 
the FRB emission detected by CHIME and GBT, respectively. 
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Extended Data Figure 5: Two-dimensional Z2-test searching result between 
two glitches We performed search in the a. 3-79 keV, b. 3-5 keV, c. 5-10 keV, 
d. 10-20 keV, and e. 20-79 keV bands using NuSTAR data. X- and y-axes 
of panels b-d are the same as that of panel a. The color map of each panel 
denotes the Z2 values, while the peak is marked with the green plus sign. The 
green contours denote Z (max) — 2.3, Z2(max) — 6, and Z2(max) — 15 where 
Z2(max) is the maximum Z3 value of the peak. 
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Extended Data Figure 6: Spectral fitting results for burst emission from 
SGR, 1935+2154 at the epochs between two glitches (Epochs B, C, and D). 
The data are extracted from NICER (black), NuSTAR FPMA(red), and NuS- 
TAR FPMB (green). All the detected burst events are accumulated. The total 
exposure of the NuSTAR data is 824 s. The accumulated spectrum is fitted 
by a model with an absorbed blackbody plus a power law with an exponential 
roll-off. Top, middle, and bottom panels show the count spectral components, 
fitting residuals, and spectra in vF, form, respectively. 


